Analisi dei dati effettuata con la suite CONTINUUM ANALITICS https://www.continuum.io/
basata sui seguenti moduli python:
tutti i software utilizzati sono open source
In [41]:
%matplotlib inline
#importo le librerie
import pandas as pd
import os
from __future__ import print_function,division
import numpy as np
import seaborn as sns
os.environ["NLS_LANG"] = "ITALIAN_ITALY.UTF8"
In [42]:
#importo il file con i dati
path=r"D:\d\05 Lavscien\autoimmunita\corr_thibya\compar_thibya_brahms.csv"
database=pd.read_csv(path,sep=';',usecols=[1, 2, 3,4,5])#colonne da utilizzare
database['valore_cap']=database['valore_cap'].apply(lambda x: round(x,2))
database.drop(['codificato','accettazione'],axis=1,inplace=True)
database.tail(6)
Out[42]:
Varibili d'ambiete in comune
In [43]:
#variabili d'ambiente comuni
cutoff_cap=2.9 #tre 2.9 r 3.3 dubbi
#cutoff_cap=3.3
cutoff_rout=1 #brahms 1-1.5 dubbi
METODO_ROUTINE="Brahms Trak Human LIA"
#METODO_ROUTINE="Siemens Immulite 2000 Chemil."
CAP="Thermo Fisher ELIA anti-TSH-R Cap250 "
In [44]:
database['cap_PN']=(database['valore_cap']>=cutoff_cap)
database['rut_PN']=(database['valore_rut']>=cutoff_rout)
database['cap_PN'].replace([True,False],['Pos','Neg'],inplace=True)
database['rut_PN'].replace([True,False],['Pos','Neg'],inplace=True)
database.head()
Out[44]:
In [45]:
database.describe()
Out[45]:
In [46]:
#sci.py moduli
from scipy.stats import chi2_contingency, fisher_exact
pd.crosstab(database.cap_PN,database.rut_PN)
Out[46]:
In [47]:
ax=pd.crosstab(database.cap_PN,database.rut_PN).plot(kind='bar',stacked=True, )
ax.legend(['Neg','Pos'])
ax.set_xlabel(CAP)
Out[47]:
In [48]:
# test chi square
chi2, pvalue, dof, ex = chi2_contingency(pd.crosstab(database.cap_PN,database.rut_PN))
print ('valore di p:{}'.format(pvalue))
In [49]:
# test esatto di Fisher
oddsratio, pvalue =fisher_exact(pd.crosstab(database.cap_PN,database.rut_PN))
print ('valore di p:{}'.format(pvalue))
In [50]:
from statsmodels.sandbox.stats.runs import mcnemar
stat,p=mcnemar(pd.crosstab(database.cap_PN,database.rut_PN))
print("valore di p:{}".format(p))
-
In [51]:
# grafico di dispersione
import matplotlib.pyplot as plt
fig = plt.figure()
fig.suptitle('Scatterplot', fontsize=14, fontweight='bold')
ax = fig.add_subplot(111)
ax.set_xlabel(METODO_ROUTINE)
ax.set_ylabel(CAP)
ax.plot(database.valore_rut,database.valore_cap,'o')
plt.show()
eseguiamo ora lo studio di regressione con tre modelli diversi
Moduli statmodels e scipy
In [52]:
# con statmodel : regressione minimi quadrati
##res_ols = sm.OLS(y, statsmodels.tools.add_constant(X)).fit() per vecchia versione
import statsmodels.api as sm
#sm.OLS(Y,X)
X = sm.add_constant(database.valore_rut )
modello_minquad=sm.OLS(database.valore_cap,X)
regressione_minquad=modello_minquad.fit()
regressione_minquad.summary()
Out[52]:
In [53]:
# con statmodel : regressione robusta (Robust Linear Model)
X = sm.add_constant(database.valore_rut)
modello=sm.RLM(database.valore_cap,X)
regressione_robusta=modello.fit()
regressione_robusta.summary()
Out[53]:
In [54]:
#importo la librearia seborn per una migliore visualizzazione grafica
sns.set(color_codes=True)
ax = sns.regplot(x=database.valore_rut,y=database.valore_cap, color="g",robust=True)
ax = sns.regplot(x=database.valore_rut,y=database.valore_cap, color="b")
ax.set_title('Regressione lineare OLS + RLM ')
ax.set_xlabel(METODO_ROUTINE)
ax.set_ylabel(CAP)
ax.set(ylim=(0, None))
ax.set(xlim=(0, None))
Out[54]:
In [55]:
sns.set(color_codes=True)
ax2 = sns.regplot(x=database.valore_rut,y=database.valore_cap, color="g",robust=True)
ax2 = sns.regplot(x=database.valore_rut,y=database.valore_cap, color="b")
ax2.set_title('Regressione lineare OLS + RLM ')
ax2.set_xlabel(METODO_ROUTINE)
ax2.set_ylabel(CAP)
ax2.set(ylim=(0, 20))
ax2.set(xlim=(0, 8))
Out[55]:
In [56]:
ax=sns.jointplot(x=database.valore_rut,y=database.valore_cap, kind="reg");
ax.set_axis_labels(METODO_ROUTINE,CAP)
Out[56]:
Ortogonal Distance Regression (Deming Regression)
In [57]:
# regressione ODR (ortogonal distance regression Deming)
import scipy.odr as odr
#modello di fitting
def funzione(B,x):
return B[0]*x+B[1]
linear= odr.Model(funzione)
variabili=odr.Data(database.valore_rut,database.valore_cap)
regressione_ortogonale=odr.ODR(variabili,linear,beta0=[1., 2.])
output=regressione_ortogonale.run()
#print (odr.Model)
output.pprint()
In [58]:
print("coefficente angolare: {ang}, Intercetta: {int}".format(ang=output.beta[0],int=output.beta[1]))
In [59]:
database_b=database
database_b['bias']=database['valore_rut']-database['valore_cap']
database_b.head(5)
Out[59]:
In [60]:
sns.distplot(database_b.bias)
Out[60]:
In [61]:
database.describe()
Out[61]:
Creo colonne con Positivo negativo dubbio in base ai cut off secificati dalle ditte
In [62]:
def discret_cap(x):
if x<2.9:
return 'N'
elif x>=3.3:
return 'P'
else:
return 'D'
def discret_bra(x):
if x<1:
return 'N'
elif x>=1.5:
return 'P'
else:
return 'D'
database['cap_PND']=database['valore_cap'].apply(discret_cap)
database['rut_PND']=database['valore_rut'].apply(discret_bra)
In [63]:
database.head(12)
Out[63]:
In [64]:
pd.crosstab(database.cap_PND,database.rut_PND)
Out[64]:
In [65]:
ax=pd.crosstab(database.cap_PND,database.rut_PND).plot(kind='bar',stacked=True, )
ax.legend(['Bra Dub','Bra Neg','Bra Pos'])
ax.set_xlabel(CAP)
Out[65]:
Creo una colonna che assume valore positivo solo nel caso in cui i due metodi abbiano dato valore opposto N con P o vieceversa)
In [66]:
def no_match(x):
if (x['cap_PND']==x['rut_PND']or x['cap_PND']=='D' or x['rut_PND']=='D'):
return 0
else:
return 1
#df.apply(lambda row: my_test(row['a'], row['c']), axis=1)
database['mismatch']=database.apply(no_match,axis=1)
#database['valore_cap'].apply(discret_cap)
per_disc=round(100*database['mismatch'].sum()/database['mismatch'].count(),2)
database.head(20)
Out[66]:
In [73]:
print("classificazioni deiverse: {} su un totale di {}".format(database['mismatch'].sum(),database['mismatch'].count()))
print(" pecentuale di discordanza e dello {}%".format(per_disc))
In [ ]: